Relaxation of the Ising spin system coupled to a bosonic bath and the time dependent mean field equation

The Ising model does not have strictly defined dynamics, only a spectrum. There are different ways to equip it with time dependence, e.g., the Glauber or the Kawasaki dynamics, which are both stochastic, but it means there is a master equation that can also describe their dynamics. These equations can be derived from the Redfield equation, where the spin system is weakly coupled to a bosonic bath. In this paper, we investigate the temperature dependence of the relaxation time of a Glauber-type master equation, especially in the case of the fully connected, uniform Ising model. The finite-size effects were analyzed with a reduced master equation and the thermodynamic limit with a time-dependent mean field equation.


Introduction
Spin models are versatile because they are simple, yet able to demonstrate fundamental phenomena, like phase transition [1][2][3]. The spin models originated from solid state physics, where the interaction between the electron spins can be due to direct exchange [4], indirect exchange [5], superexchange [6] or double exchange [7], but many other complex physical systems can be modelled using a simple Ising or Heisenberg model, like nuclear spins [8,9], and even social situations [10]. It is also important in modern applied physics since one branch of adiabatic quantum computers-e.g. the D-Wave system [11]-are based on finding the global minimum of an artificial spin system [12,13].
The Ising model is defined via its energy or, in the quantum case, where it is called Heisenberg model via its Hamiltonian operator. All the equilibrium properties can be determined from these quantities, e.g. heat capacity, magnetic moment, susceptibility etc., since they can be derived from the partition function. On the other hand to calculate inequilibrium properties like the relaxation time, we need to know the dynamical equation of the system. The Ising model does not have a natural dynamics, and although to the Heisenberg model we can associate the Schrödinger or the Heisenberg equation, it will only generate a unitary time evolution. Therefore it cannot be responsible for a final, thermal distribution. For this we need some interaction with the surrounding environment. Close to thermal equilibrium it is often assumed, that the dynamics is driven by the gradient of a mean field free energy, since at equilibrium it must be zero [14][15][16][17]. This phenomenological approach can account for some dynamical critical behaviour, like the dynamical slowing down, but like every mean field approximation it is valid only if the fluctuations are small. For example in one dimension a master equation gives a better description [18].
A more fundamental approach is to derive the effective dynamical equations from a system-plus-bath model, where the bath is a thermal reservoir. This method let us track how the parameters of the macroscopic equation depends on microscopic properties. To describe such a system we must use the tools of open quantum systems [19,20] like the Redfield [21] and the Lindblad equation [22]. These equations have countless applications in quantum biology [23,24], quantum optics [19], cold atomic gases [25], chemical physics [26] as well as it also being relevant in quantum computing [27][28][29].
In solid matter the interaction with the phonons is always present and the electrons, as charged particles are coupled to the photons, thus we can assume that the spin system is in a bosonic bath. Quantum dissipation and relaxation of spin systems in a bosonic bath as well as in magnetic field have been investigated by many authors. [30][31][32][33].
In this paper the interaction between an adiabatic computer and its environment is meant to be small, so the weak coupling Lindblad equation is be used, but of course, there are improved methods to describe open quantum systems, like slippage initial condition [34,35], the Nakajima-Zwanzig equation [36,37] or the polaron transformation [38].
The structure of this paper is the following. In section 2 we briefly summarize the derivation of a Glauber-type master equation [39] based on the Redfield equation. Using this microscopic approach we can see how the Fourier transform of the bath correlation function appears in the final master equation. In section 3 we investigate the temperature dependence of the eigenvalues of the transition matrix because they contain relevant information on the time scales of the system, e.g., the relaxation time. We also give an upper bound to the smallest nonzero eigenvalue. In section 4 the dynamics of the uniform, fully connected Ising model is investigated, and we show that the relaxation time diverges in the thermodynamic limit as the temperature approaches the critical temperature. In section 5 a time dependent mean field equation is revisited, which is used in section 6 to extend the previous results to infinite sizes.

Master equation of quantum Ising system
In general if a system is connected to a bath, then its Hamiltonian operator is where H acts only on the system of interest, H B only on the bath, and H I is the interaction between the two subsystems, it can be written as where A α and B α are system and bath operators respectively. The dynamics of the total system is described by the von Neumann equation.
After the Born and the Markov approximation an effective equation can be derived for the density operator of the system of interest. (r≔Tr Br tot ). where is in interaction picture, C ab ðtÞ ¼ hB I a ðtÞB b i B is the bath correlation function, and h.c means hermitian conjugate. This is the Redfield equation in weak-coupling limit [21]. The Born approximation is valid if the interaction between the system and the bath is small, and we can use the Markov approximation if the relaxation time of the system of interest is much larger than the decaying time of the bath correlation function. After the so called secular or rotating wave approximation [40] one can get to the Lindblad equation [19,22]: where A a ðoÞ ¼ P ij jiihijA a jjihjjd o;ε j À ε i and |ii is the eigenvector of H with eigenvalue ε i . H LS is the Lamb shift Hamiltonian, and γ αβ (ω) is the Fourier transform of the bath correlation function.
There are two standard bosonic baths: the bath of phonons and the bath of photons. For pho- 1À e À bo , which is called Ohmic case and for photons g sup ðoÞ � o 3 1À e À bo , which is a super-Ohmic case. The frequency ω c is the cutoff frequency. If we assume, that ω c is large compared to the energy distances of the system, then e À o jo c j � 1. Fig 1 shows the main features of the two γ functions. The main difference is that γ ohm is strictly increasing and γ ohm (ω = 0) * k B T, but γ sup is non-monotonic, and γ sup (0) = 0. In the easiest case γ αβ / δ αβ .
The advantage of the weak-coupling limit is that a master equation can be derived for the diagonal elements of ρ. where The system converges to the Boltzmann distribution if W ij satisfies the detailed balance condition i.e., W ij = W ji exp(−β(ε i − ε j )). Both the Ohmic and the super-Ohmic bath satisfy it, because gðÀ oÞ ¼ e À bo gðoÞ ð7Þ If the system of interest is the Ising model, then a master equation can be derived [41], which contains only the diagonal elements of the density matrix. The Hamiltonian is where s z i is the Pauli z-matrix and the corresponding eigenvectors are with eigenenergies The easiest way to couple the system to the bath is via a Pauli matrix i.e., A a ! s x i . Using σ z in the interaction instead of σ x would not give any relevant dynamics since the system and the interaction Hamiltonians would commute. The peculiarity of this system is that the populations decouple even without the secular approximation.
The s x i operator acting on |Si only flips the ith spin, so the W SS 0 matrix element is With Eqs (6) and (11) we have a dynamics for the Ising model.
This matrix is temperature dependent, and it has at least one zero eigenvalue, which is the eigenvalue of the equilibrium distribution: For constant temperature the general solution of (12) is where P R m s are the right, and P L m s are the left eigenvectors of M with −λ μ eigenvalues. All the λ μ s are nonnegative. If the system is ergodic, then there is only one zero eigenvalue, and the other λs are positive. Let the smallest positive be λ min and the largest be λ max . The relaxation time is t r = 1/λ min . This is the time scale in which all but the equilibrium mode dies out. The other relevant time scale is 1/λ max , which is the characteristic time of the fastest mode. If, for example, this spin system is a quantum computer, then the fastest mode is more important because if the computation is slower than this time scale, then the environment is not negligible. In other words, λ min is important if we want the system to relax thermally, and λ max is important if we want to avoid any thermal influence.

Temperature dependence of the eigenvalues
Both the smallest and the largest eigenvalues carry relevant information, and since M(β) is temperature dependent λ min (β) and λ max (β) are too.
At high temperatures we can determine the temperature dependence of all λs by simply Taylor expanding γ(ω; β) for small β.
where α = 1 in the Ohmic, and α = 3 in the super-Ohmic case. The transition matrix inherits this temperature dependence: M SS 0 (β) * β −1 , and hence λ * β −1 . Despite the high-temperature limit, where the elements of the dynamical matrix M diverges, in the low-temperature limit they converge.
It means all the eigenvalues also converge. As a consequence, we cannot slow down arbitrary all the modes by reducing the temperature. We have an upper limit in time for quantum computing. Of course, this calculation is valid only for a time-independent system, but the main features apply to more general cases. Without an external magnetic field (h = 0) at zero temperature, the equilibrium Boltzmann distribution prefers only the two spin configurations with the lowest energies: where S g and −S g are the ground states. At zero temperature, there is one more eigenvector with a zero eigenvalue: The question is how λ min (β) behaves at low temperature. We can give an upper bound. First let us introduce the following symmetric matrix: This transformation doesn't affect the eigenvalues, and the eigenvectors transform likẽ SinceM is symmetric its right and left eigenvectors are the same, and now the variational method applies to it: whereP is an arbitrary vector with P SP 2 S ¼ 1, and it must be perpendicular to the equilibrium vector (P eq S � ffi ffi ffi ffi ffi ffi P eq S q ), because λ min is the second smallest eigenvalue ofM. Let Then according to (21) l min � À 1 2 ðM S g ;S g ÀM S g ;À S g ÀM À S g ;S g þM À S g ;À S g Þ where d(S, S g ) is the Hamming distance, and the S 7 ! −S symmetry was used. In the bosonic bath where DE S ≔E S À E Sg > 0. At low temperatures this is the sum of some e À bDE S functions, so λ min (β) can be estimated from above with an exponential function. � 0:44). The left figure is in log-log scale, where we can see, that at low temperature the eigenvalues has a β −1 temperature dependence, and λ max (β) converges, and the right figure with lin-log scale shows, that λ min (β) goes to zero exponentially.

Eigenvalues of the uniform Ising model
The M matrices are 2 N × 2 N large; therefore, we cannot see how the eigenvalues behave at the thermodynamic limit. However, the uniform, fully connected Ising model is so symmetric that an effective equation can be derived with the same relaxation time as the original equation.
The energy of the model is The 1/N factor is to keep the energy extensive and J > 0. Given an S microstate, it consists of N " spins with S i = 1 and N # spins with S i = −1. The number of spins is constant, i.e., N " + N # = N = fix. The energy of such a configuration is If N is fixed, then the energy is the function of only N " . The symmetry of the system is that we can permute the spins in any way, the energy and the M matrix remains the same. If in the dynamics the initial condition also has this symmetry, then the P S will inherit this property. The slowest mode propagates between the two deepest valleys of the energy landscape, which are the " " . . ." and ## . . .#. Assume that initially P ##. . .# (t = 0) = 1, and we want to determine relaxation time, where P ##. . .# (t r ) � P " ". . ." (t r ). Since both the equations and the initial condition have the permutation symmetry all the probabilities, which have the same up spin have the same value, e.g. for 3 spins P "## (t) = P #"# (t) = P ##" (t) 8t. The probability can only flow between spin configurations if the Hamming distance between them is 1. Let us introduce the following probabilities: where the prime denotes that only such configurations count where there are N " up spin. We can give a closed set of differential equations which only contain these new P N " probabilities.
This master equation has only N + 1 variables instead of 2 N , thus is easy to simulate for large systems. A comparison between the quantum and the thermal simulated annealing of the fully connected Ising model was investigated by Wauters et al. using a similarly reduced master equation [42]. Eq (27) has the form and we want to determine the lowest (nonzero) eigenvalue of M red , which is the same as the lowest (nonzero) eigenvalue of M. The matrix M red is sparse, because it is a tridiagonal matrix, i.e., only the main diagonal, the first diagonal below and above the main diagonal are nonzero. Fig 3a shows the temperature dependence of λ min for different system sizes. As N increases we can see, that around the critical temperature (which is β c J = 1) the behaviour of the system changes. At Fig 3b we can see better that above the critical temperature (T > T c ) for large N values λ min converges, meaning for every system size there is a finite relaxation time. At the critical temperature (T = T c ), it follows a power law (λ min / N −0.5 ). Below the critical temperature (T < T c ) λ min goes to 0 for large N, but doesn't follow a power law. This behaviour is the famous critical slowing down phenomenon. From the N ! 1 thermodynamic limit we can determine the dynamical critical exponent. Fig 4 shows λ min (T, N ! 1) as a function of the reduced temperature ( TÀ T c T c ). This follows an easy power law, because λ min / T − T c . In the next section we will see that this result can be obtained from the mean field approximation.

Time dependent mean field equation
Since the primary interest is the magnetization (m i ≔ hS i i), we would like to derive a differential equation for it. Penrose showed that with some approximations, this is possible [43], and as the number of neighbors increases, this becomes more and more accurate.
Using the definition of m i and the master equation we get The W S 0 S matrix component is nonzero if the Hamming distance between S 0 and S is one. Introducing we can rewrite the double sum in (30).
In the second step the (Λ i (S, n) − S i ) = 2S i δ in identity was used. The nonzero elements of W are a function of the energy difference: whereh i ¼ P j J ij S j þ h i , so this is still the function of the S random variable, but because J ii = 0 it is not a function of S i . Since S i can be only 1 or −1 the gðÀ 2h i S i Þ as a function of S i must have the form. Using (7) yields then substituting back to (32) gives Eq (36) is similar to the Callen equation [44,45]ðhS i i ¼ htanhðbh i ÞiÞ, where the averaging is outside the hyperbolic function. In order to get a closed equation for the expected values, the average must move inside, and instead of the S i random variables, their m i expected values must be written.
This mean field approximation is valid as long as the fluctuations are small, which holds if every spin interacts with many other spins. This can be due to long range interaction or in high spatial dimensions [46]. The right-hand side of Eq (37) contains the self-consistent equation from the equilibrium statistical physics; hence if the equation of state is satisfied, then _ m i ¼ 0. Eq (37) contains both the real-time and the temperature of the bath. The temperature can also be time-dependent, and in that case, we could get thermal annealing, but if the temperature is constant, we can determine the relaxation time and the dynamical critical exponent. If m(t) = m eq + δm(t), where m eq is the equilibrium solution and δm(t) is small, then the linearized equation of (37) is Using the 1 cosh 2 ðxÞ � 1 À tanh 2 ðxÞ identity, and the equation of state we get to Eq (40) contains the inverse susceptibility of the mean field Ising model.
Substituting the inverse susceptibility back into (40) yields This is a well-known equation in the theory of dynamical critical phenomena [47], but it is usually derived from the _ m ¼ À G@ m F MFA phenomenological equation. Now we can see, how it is related to a master equation and the spin-boson model. If the system is symmetric in a sense, that all the spins behave the same, then (43) simplifies to where w À 1 ¼ P j w À 1 ij .

Time dependent mean field equation for the uniform Ising model
As before in section 4 the uniform Ising model will be studied because in the equilibrium case in the thermodynamic limit, it gives back the exact results. According to (42) the mean field free energy is and the time dependent mean field equation is If h = 0 the critical temperature is T c = J, and above this temperature the equilibrium solution is m eq = 0. The inverse susceptibility is therefore where Γ = 2γ ohm (0; β)β = 2η in the Ohmic bath. In the super-Ohmic bath, this is zero because γ sup (0) = 0, i.e., there is no transition between states with the same energy. In a real physical system it means that besides the Glauber dynamics some other effects are not negligible, e.g. 2 spin flips. Eq (48) is the same result that we have already seen in Fig 4. As in equilibrium statistical physics, the mean field approximation gives back the exact result for the uniform model in the thermodynamic limit. At the critical temperature the inverse susceptibility is zero, the linear term vanishes, and we need the higher order terms. Taylor expanding (46) at T = T c around m = m eq � 0 up to the third order gives.
which has the solution for large ts, which means there isn't a characteristic time.
Below the critical temperature, the linearized equation is good again; only m eq changes. On the other hand Eq (46) can give a different solution from the (29) master equation. If we want to compare these two equations, the initial condition must also be the same, which gives a restriction to the initial condition of (29). In the mean field approximation, the probability is a product of the one-particle probabilities: which means the initial probability of (29) is If h = 0, then the master equation must converge to the m = 0 solution, but the time-dependent mean field equation finds one of the minima of the free energy if initially m 6 ¼ 0. If h is finite and m is in the valley of the global minimum (point A in Fig 5), then in the N ! 1 limit the solution of the master equation converges to the mean field solution (Fig 6).
On the other hand, if initially, m is in the valley of the local minimum (point B in Fig 5), then the solution of the mean field equation converges to this local minimum, but the solutions distribution, so they tend to approach m eq for large N values in the t ! 1 limit, but as N grows, so does the relaxation time. In the thermodynamic limit the relaxation time diverges as in Fig 3. In the N ! 1 limit the solution of the master equation converges to the solution of the mean field equation and none of them will approach the global minimum, because the relaxation time will be infinite.

Conclusion
In this work, we have studied a Glauber-type master equation derived from the spin-boson model. The most relevant dynamical properties are encoded in the eigenvalues of the transition matrix of the master equation. They are temperature dependent and behave significantly differently below, above, and at the critical temperature as a function of the system size. In the case of the uniform, fully connected Ising model, in the thermodynamic limit, above the critical temperature the relaxation time follows a power law: t r / (T − T c ) −1 .
Using the time-dependent mean field equation, we could also investigate the thermodynamic limit. Its dynamics differ significantly from the finite size master equation if the initial condition is close to a local minimum of the free energy, which means that the relaxation time diverges. If it is close to the global minimum, then the solution of the finite size master equation converges to the solution of the time-dependent mean field equation.